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LIST OF MAIN DEFINITIONS 


The notation for orbital elements is standard. Figure 2-1 
presents the nomenclature of the coordinate system. In addition, the 
following are the main definitions used: 

Contraction of Orbit 

X = p ae = ae/ H = semi-major axis X eccentricity/ scale height 
z = a/ a = ratio of semi-major axis to initial semi-major axis 

Ballistic Entry 

Modified Chapman Variables 

u = V^'cos^' y/ gr J 

V = V^/ gr = dimensionless velocity squared 

Y = 2Z 

$ = - <s/ p r' sin y 

X = - Log V 

"y = 2Z / 

$ = - sin y 

•n = - ZZ/'s/^ sin y . 

S = sin y. / sin y 

1 

— ^i 

V. = V. e 

1 1 
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where 


p = atmospheric density 

S = reference surface area of space vehicle 
= drag coefficient of space vehicle 
m = mass of space vehicle 

r = radial distance from center of planet to space vehicle 
P = 1/H = inverse scale height of atmosphere 

V = absolute velocity of space vehicle 

Y = flight path angle 

g = gravitational acceleration at radial distance r 
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ABSTRACT 


A space object traveling through an atnnosphere is governed by 
two forces: aerodynamic and gravitational. On this premise, equations 
of motion are derived to provide a set of universal entry equations 
applicable to all regimes of atmospheric flight from orbital motion 
under the dissipative force of drag through the dynamic phase of 
reentry, and finally to the point of contact with the planetary surface. 

Rigorous mathematical techniques, such as averaging, Poincare's 
method of small parameters, and Lagrange's expansion, are applied to 
obtain a highly accurate, purely analytic theory for orbit contraction and 
ballistic entry into planetary atmospheres. 

The theory has a wide range of applications to modern problems 
including orbit decay of artificial satellites, atmospheric capture of 
planetary probes, atmospheric grazing, and ballistic reentry of manned 
and unmanned space vehicles. 
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CHAPTER 1 


INTRODUCTION 

Since earliest recorded times man has always been fascinated 
by the motions of objects in the sky. This natural curiosity resulted 
in the development of astronomy as the first science which in turn 
acted as a catalyst for the advancement of mathematics. Many of 
the classical scholars such as Ptolemy, Copernicus, Kepler and 
Newton were both astronomers and mathematicians, and all of them 
made important contributions to a field now known as celestial 
mechanics. 

Although Newton foresaw the possibility of artificial satellites 
and even computed the decay due to a resisting medium, it has only 
been in recent times that man has had sufficient knowledge of the 
atmosphere to develop adequate theories of satellites in low orbits. 
With empirical data from the first artificial satellites, King-Hele 
(1964) was able to publish a monograph presenting a large step in 
sophistication of analytic theories. He, and others, were primarily 
interested in applying classical celestial mechanics techniques in 
studying the slowly varying orbital elements of the osculating orbit, 
under the perturbation of atmospheric drag. Since the original 


1 



artificial satellites were not intended for recovery, the researchers 
confined their theories to the prediction of the satellite's orbital 
motion until its eventual burn-up in the Earth's atmosphere. 

With the advent of manned space exploration and recoverable 
space probes, it became necessary to develop accurate theories of 
the entry phase, during which the altitude, velocity, deceleration, 
and the heating rate all vary rapidly. Space dynamicists such as 
Chapman, Loh and Yaroshevskii produced analytic theories to 
describe the reentry of ballistic and lifting spacecraft. In these 
analyses strong physical assumptions were made and soon the 
equations of the reentry scientists became very different from those 
of the celestial mechanics scholars. 

At the University of Michigan the theories of King-Hele, 
Chapman, Loh and Yaroshevskii have been reviewed with an aim of 
achieving greater generality and flexibility in application. Vinh 
and Brace (1974) removed Chapman's restrictive assumptions by 
introducing modified Chapman variables which can be applied to 
three dimensional entry trajectories. Since the problem of orbit 
decay and the problem of reentry deal with the same physical 
phenomenon, namely, flight of a space object in an atmosphere, 
they must obey the same equations. Vinh, Busemann and Culp 
(1975) developed the set of exact universal entry equations 
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applicable to spacecraft in all regimes of flight from orbital motion, 
entry phase, glide and touchdown. The theses of Brace (1974) and 
Bletsos (197 6) clearly demonstrate the superiority of the more 
general approaches. The purpose of this report is to illus- 
trate the application of the universal entry equations to orbit decay, 
as well as ballistic entry into a planetary atmosphere and to present 
rigorous mathematical analyses of the resulting math models. In 
order to provide a firm foundation the simplest cases were assumed 
in the math model: exponential atmosphere, constant por [3 r, and 
spherical rotating planet. However, the techniques can be extended 
to include the effects of varying scale height and oblate atmosphere. 
The report is organized into chapters as follows. 

In Chapter 2 the universal entry equations are derived from 
aerodynamic equations of motion. Next, the perturbation equations 
for satellite motion in a rotating atmosphere are deduced from the 
universal equations. Then, in Chapter 3, the variational equations 
for orbit decay are derived by applying the averaging technique to 
the perturbation equations. The solutions of these equations are 
presented in Chapter 4 as the theory of orbit decay. Poincare's 
method of small parameters is applied to the basic nonlinear first 
order differential equation for contraction of orbits of small and 
moderate eccentricity and exact solutions are obtained for the first 
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five order terms to find the semi-major axis and eccentricity in 
parametric form. Lagrange's expansion is used to obtain the semi- 
major axis as a function of eccentricity. For highly eccentric orbits, 
the governing differential equation is derived directly from the basic 
equation and solved in closed form parametrically. First order 
solutions are found for the other orbital elements, giving the orient- 
ation of the orbit. Solutions are also derived for circular and nearly 
circular orbits. For the main effect of orbit contraction, a uniformly 
valid, mathematically rigorous theory for all eccentricities, 0 < e < 1, 
is presented. 

In Chapter 5 the problem of time in orbit is considered. 
Difficulties arise in obtaining exact integrals, but an approximate 
method gives excellent numerical results. The time solution is used 
for the time at any point during the orbit decay, as well as for the 
maximum lifetime of the satellite for the assumed atmospheric model. 

In Chapter 6 a higher order theory for ballistic entry is 
developed by applying Poincare's method to a system of nonlinear 
differential equations. The problem is divided into three separate 
cases: entry from circular orbit, entry at small initial flight path 
angles and entry at moderate and large flight path angles. The second 
order solution of the first case gives a significant improvement over 
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the Yaroshevskii solution which enters as the zero order term. 
Similarly, the solution to the third case shows marked improvement 
over the classical solution of Chapman, again entering as the zero 
order term. The second case bridges the gap between the first and 
third cases and provides adequate numerical results. 

In Chapter 7 conclusions are drawn and recommendations 
for future research are made. 
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CHAPTER 2 


THE UNIVERSAL EQUATIONS FOR ORBIT DECAY 

AND REENTRY 


2. 1 Forces on a Satellite in Orbit 

Only two forces will be considered to effect the motion of a 
space object in an atmosphere: gravitational and aerodynamic. It 
is assumed that the planet is spherical with uniform mass distribu- 
tion, so that the gravitational force is an inverse square force of 
attraction with acceleration 


g(r) = 


( 2 - 1 ) 


where r is the radial distance from the center of the planet and p, is 
a positive constant which, for two-body relative motion is G(M + m) 
where G is Newton's universal gravitational constant, M is the mass 
of the planet and m is the mass of the spacecraft. The atmospheric 
force is in the form of drag, D , which acts in a direction opposite 
to the velocity relative to the ambient air 


^A = IP^S^A 


(2-2) 


where is the drag coefficient for a reference area S and p is the 
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atmospheric density. The atmospheric density is assumed to vary 
exponentially with altitude 


P(r -r) 

Po 

p = p e 

Po 


where p is the density at the initial periapsis, r 

P P 

o o 


constant 


P 


_1 

H 


(2-3) 


and j3 is a 


(2-4) 


where H is the scale height. For simplicity, it is assumed that 
there is no lift force. 

These forces represent perhaps the most elementary of 
realistic models. The emphasis will be to provide a pure mathemat- 
ical treatment of the orbit decay problem which will serve as a firm 
foundation for the development of more sophisticated solutions which 
include oblateness of the atmosphere and variation of scale height, 
as well as other important features. 

The equations of motion are written with respect to an inertial 
reference frame with the origin at the center of the planet. In Fig. 
2-1, V is the absolute velocity of the satellite 

V = (2-5) 

where V is the velocity at the point M, of the ambient air relative 
to the planet center. If w is the angular velocity of the rotating 
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atmosphere, then 


V = rw cos d> (2-6) 

e 

where 4> is the latitude of the point M. 

If 4^ ' is the angle between and V, then by the cosine rule 
2 2 2 ♦ 

V. = V + V -2VV cosij;' . (2-7) 

A e e 

The heading angle, , is the angle between V and the projection 
V of V on the horizontal plane and is related to the latitude (j) and 
the inclination i of the orbital plane by the well-known relation 


cos 4^ cos 4> = cos i . (2-8) 

The flight path angle, y , is the angle between V and V and is 

H 

related to 4< s-nd 4^ ' by 

cos 4 j ’ = cos 4 j cos y (2. 9) 


Near periapsis , where the aerodynamic drag is most effective, the 
satellite travels in a nearly horizontal direction and hence y is 
small, allowing the approximation 


V cos ijj ' = V cos4j cosy « rwcoscj) cos4j = rwcos i . (2-10) 


Substituting Eqs. (2-6) and (2-10) into Eq. (2-7) gives 

2 2 

2 2 2rw . r w 

V " ^ ^ ^ — T" • 


( 2 - 11 ) 


The rotation of the atmosphere is so slow that the term w can be 
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neglected. For the small term rw/ V, it is appropriate to use an 

average value. King-Hele suggested using the average value at 

perigee r / V to replace r/ V. Since the inclination i usually 
o o 

varies by less than 0. 3 during a satellite's lifetime, it may be 

taken equal to its initial value i . The result is King-Hele 's 

o 

expression (King-Hele, 1964) 


v/ = £ 

A 


( 2 - 12 ) 


with the slightly different average constant 

2r w 

f = ( 1 - — cos i ) 

V o 


(2-13) 


Thus, in terms of the absolute speed, the drag force is 

which acts opposite to the direction of the velocity V of the 
satellite relative to the ambient air. 


2. 2 The Universal Equations of Motion 


For the flight of an aerodynamic vehicle in a nonrotating 

atmosphere with a lift coefficient and a drag coefficient , 

L ® D 

it is customary to use the equations of motion with the notation 
of Fig. 2- 1 . 
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V sin -y 

V cos y cos 4^ 
r cos <j) 

V cos y siniji 


psc^v 


( 2 - 15 ) 


g sinv 


p SC^V cos O' yZ 

^ — - ( g - — ) cos V 


p SC V sinr 

Xj 

2m cos Y 


cos Y cos i[j tanej) 


where the bank angle cr is defined as the angle between the local 
vertical plane containing the velocity and the plane containing the 
velocity and the aerodynamic force. 


Using the modified Chapman variables 


V cos Y 


,Z = 




( 2 - 16 ) 


and the dimensionless independent variable 


r V 

s = J ^ y 


( 2 - 17 ) 


the exact universal equations for entry trajectories into a planetary 
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atmosphere assumed to be at rest are derived (Vinh et al. , 1975) 


dZ 

ds 

du 

ds 


lx 

ds 


li 

ds 


dc)) 


ds 


ds 


- (3 r ( - — + -i- ) Z tan V 

PP dr Zpr 2(3^ dr 


zyfTr 


Pr Zu , , , L 

{ 1 + — coscr tanv + 

COSY ‘ C * 


sin Y 


VPr Z 
COS Y 

cos 
cos (j) 


D 


L cos Y 

“ — coscr + *- 


D 




2^pr Z 
( 1 - 


2 

CO s Y 


r Z 


u 


(2-18) 


sin i)j 

VTr 2 
2 

cos Y 


( — ^ sin cr 


cos 




Y cos i|j tancj} ^ 


2. 3 Satellite Motion in a Rotating Atmosphere 

Next, the universal entry equations will be transformed to 
obtain the equations for satellite motion inside a rotating atmosphere. 
Using the following relations from spherical trigonometry 
cos cj> cos cjj = cos i 

cos ^ sin i|j = sin i cos a (2-19) 

cos a = cos (|) cos ( 0 - S2) 

the last three equations of (2-18) are transformed to 
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do- 

ds 

dn 

ds 


d i 
ds 


1 - 


VTr 


Z sin a , L , . 

( — ) sin 0 - 


tan i cos -y 

Z sin a 

2 C 

sin i cos Y D 

Vp r Z cos a ^L. 

cos Y D 


D 


• ) sin (T 


( 2 - 20 ) 


) sin cr 


where is the longitude of the ascending node of the osculating plane 
and a is the angle between the ascending node and the position vector. 

For satellite motion in a rotating atmosphere there is a 
complication in that the drag force is modified by the factor f and is 
directed opposite to the velocity V , not the absolute velocity V . 

At the same time there is the simplification of no lift. 

In Fig. (2-2) the aerodynamic force diagram used in deriving 
the equations of motion (2-15) is illustrated. In addition the velocity 

V with respect to the ambient air, and the drag force D , 

opposite in direction to are shown. For the present derivation 

the lift force L is removed and the vector drag D is replaced by 

D . The force D can be decomposed into one component in the 
orbital plane and one component normal to the orbital plane. Since 

V is small, V is nearly aligned to V and the drag component in 

© A. 

the orbital plane can be considered to be directly opposite to V, 


with magnitude D as given by Eq. (2-14). The normal component 
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of the drag, , is found from the projection of 


Sa = -iPSfc^v^ 


V 


V 


( 2 - 21 ) 


in the direction orthogonal to the orbital plane. Since V is in the 

orbital plane and V = V + V , the projection of V on the normal 

A e A 

to the orbital plane is the same as the projection of which has 
magnitude 


V sin ijj = rw cos cf) sin 4 / 
e 


= rw sin 1 cos a 


( 2 - 22 ) 


Thus, the vector has magnitude 


D 


N 


1 V 

— p Sf rw sin icosa — 


(2-23) 


From King-Hele's expression (2-12) it can be written as 


P SfCj^V 


N 


2f 


T7T 


r w sin 1 cos or 


(2-24) 


and its direction is opposite to the vector L sin o' in Fig. (2-2). 

The final result of this analysis is that, in Eqs. (2-18) and 
( 2 - 20 ) the is replaced by the modified drag coefficient 
the component cos o' is deleted and the component C_ sino' is 
replaced by 

1/2 rw 

C sino' = - f C ( ) sini cos a . (2-25) 

i-< D V 
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The modified Chapman variable Z is most effective in analyz- 


ing the entry phase of the vehicle. While the vehicle is still in 


orbit, it is used in the following form with constant p 




P (r - r) 
^o 


(2-26) 


where the dimensionless constant Z is 

o 


Z 

o 


P SfC^r 

Po ° Pn 
o o 

2m 


(2-27) 


Eqs. (18) and (19) are now rewritten introducing the equation for 
r/r to replace the equation for Z. 


_d 

ds 



du 

ds 



) tan y 


2 Z u 

, o 

- u tany - 

cos y 



r) 


ds 


2 

cos y 
u 


dg 

ds 


r wZ , p(r - r) 

p o 5/2 . . P 

, o , r , cosisingcosg o 

1 + ( — ) e 




u 


cos y 
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dn 

ds 


5/2 


d i 
ds 


r w Z 
P o 
o 





sin a CO s a 

1/2 

u cos -y 


e 


5/2 . . 2 

sinicos a 


1/2 

u cos -y 




(2-28) 


Eqs. (2-28) bridge the gap between the satellite theory and 
the entry theory. As a matter of fact, they can be used to follow 
the motion of a vehicle subject to gravitational force and drag force 
of a rotating planet for its entire life in orbit until the end of its 
entry and contact with the planetary surface. The accuracy depends 
on the readjustment, for each layer of the atmosphere, of the 
"constant value" p . These equations are most useful for analyzing 
the last few revolutions and the entry phase. In Chapter 3 the entry 
variables r, u and y will be transformed into orbital elements and 
the variational equations of orbit decay will be derived. 
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Fig. 2-1. Nomenclature. 

The space vehicle M is loca^d by the position vector 
r and has absolute velocity V. Vpj is the horizontal 
component of Vprojected through the flight path angle 
•y while and V are the velocity with respect to the 
ambient air and tSe absolute velocity of the atmosphere 
respectively. Vpj and V are in the local horizontal 
plane and form an angle 5 which is called the heading 
angle. 0 and <{) are the longitude and latitude. , the 
longitude of the ascending node, and i, the inclination, 
are orbital elements of the osculating orbit, w is the 
angular velocity of the atmosphere. 
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Fig. 2-2. Aerodynamic Forces. 

In the derivation of the equations of motion 
for a rotating atmosphere, a lifting body with 
lift L , bank angle c and drag D is compared 
with a nonlifting body with out of plane drag D . 
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CHAPTER 3 


THE VARIATIONAL EQUATIONS FOR ORBIT CONTRACTION 


3. 1 The Osculating Orbit 


The osculating orbit is the trajectory the vehicle would follow 
if the perturbing force (the drag) suddenly vanished. By setting 
= 0 in Eqs. (2-28), the equations for the osculating orbit are 
found. 


_d 

ds 




tan y 


du 

ds 


- u tan y 


±L 

ds 


2 

cos y 
u 


do; 

ds 


= 1 


dn 

ds 

di 

ds 


0 

0 


(3-1) 


Dividing the first of Eqs. (3-1) by the second and integrating 


r 


u 


(3-2) 
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By taking the derivative of the second equation in (3-1) , 


d^u 

ds^ 


+ u = 1 


which has the solution 


(3-3) 


u = 1 + cos (s - ) 


(3-4) 


Substituting Eq. (3-4) back into the - 7 — equation, it is easy to show 

ds 


that 


cos y = 


u 


2u + C^ -1 


(3-5) 


Changing the indices of the constants, the general integration is as 
follows: 


2 

cos y 


u 


2u - C 
u 


1 


u 


1 + yjl - C J COS ( S - C^) 
O' + C . 


(3-6) 


Q 


Actually there are only 5 constants since s is equivalent to a. The 
sixth constant of integration is found from the time equation (2-17). 
The first three constants of integration can be evaluated by taking 
the origin of time at the instant of passage through the periapsis. 
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2 

cos y 


2 
u 

2u - (1 - e^) 
u = 1 + ecos (o' - cj) (3-7) 

a(l - ) 

r = : r- 

1 + e cos ( O' - 00 ) 

These transformation equations provide the relation of the entry 
variables r,u and y and the orbital elements a, the semi-major axis, 
e , the eccentricity, and oo , the argument of the periapsis. 

3. 2 The Variational Equations 

While the vehicle is in orbit, is small, acting as a 
perturbation, and so the orbital elements actually vary slowly. 

Thus, the variation equations, or perturbation equations, can be 
obtained from Eqs. (3-7) by taking the derivatives and assuming that 
a, e and oo are varying quantities. First, the derivatives are taken 
with respect to s to obtain; 


de 

2Z u^ . 2 

o 1 cos y 


P (r -r) 

Po 

ds 

3 1 u 

e cos y 

) \ r / 

Po 

6 

da 

- 2 Z au^ 

o / r \ 


• r) 

ds 

.,2 3 \ r ) 

(1 - e ) cos y p^ 

0 
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da 

ds 


doj 

ds 


u tan Y 

ri 

yje - (u- 1) 


2Z u / 

rr I 

cos Y /e - (u - 1) ( 


U(u- 1) 

2 2 
e cos Y 



(3 (r - r) 

Po 

e 


(3-8) 


Along with the last three equations of (2-28) this provides all the 
derivatives with respect to s. It is interesting to note from the first 
of Eqs. (3-8) that the eccentricity does not decrease continuously 
under the influence of atmospheric drag, but while decreasing 
secularly, it oscillates between maximum and minimum values as 
the flight path angle passes through minimum and maximum values 
respectively during each orbital period. 

Finally, the variational equations will be put in terms of the 
more familiar eccentric anomaly E as the independent variable. 
Since 


and 


r = a(l-e cos E) 


dr 

dE 


ae sin E 


(3-9) 


(3-10) 


The differential relation between s and E is 


ds 

dE 


dr ds 
dE dr 


e sin E 


i 


(1 - e cos E) tan Y 


1 - e 


1 - e cos E 


(3-11) 
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From Eqs. (2-28), (3-8) and (3-11) the variational equations are 


obtained. 


/ 


de 

dE 


- 2Z (1-e^) 
o 



„ / 1 + e cos E \ 

cosE |— — — I 

» 1 - e cos E / 


1/2 ■ 

o 


r) 


2 3/2 

^ - 7 7 (1 + ecosE) ^ Po 

dE ■ o r ^.1/2 ® 

p ( 1 - e COS E) 

o 




da d(o 


dE dE 1 - e cos E 


1 + I ■ ■ I sinE (1 -e cos E) 

P„ 


1/2 


dn 

dE 


r w Z 

Po ° 




P(r -r) 

1/7 

X ( 1 + e cos E) e 


5/2 5/2 

— I (1-ecosE) 


1/2 


P (r - r) 


X (1 + ecosE) sin O’ cos or e 


o 


di 

dE 


r w Z 

Po ° 


5/2 


5/2 


V‘- 


- (7^) 

• rv * 


- e cos E) 


p(r - r) 

1/2 2 Po 

X (1 + ecosE) sin i cos a e 

(3-12) 

The variational equations apply to the space vehicle for the entire 
lifetime in orbit from the initial eccentric orbit until circularization 


and the onset of reentry. In Chapter 4 mathematically rigorous 
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analytic techniques will be employed in order to solve the problem of 
of the orbit decay phase. 
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CHAPTER 4 


THE THEORY OF ORBIT DECAY 


4. 1 The Average Equations for Orbit Decay 


The variational equations (3-12) exhibit both oscillatory or 
periodic behavior and slowly varying secular behavior. The focus 
of this chapter is to solve for the slowly varying orbital elements 
of the osculating orbit and to avoid the periodic terms. Fortunately 
there is a rigorous method, the averaging technique (Bogoliubov and 
Mitropolsky, 1961), which allows the periodic terms to be replaced 
by the average value and which converges to the correct analytic 
solution as the independent variable goes to infinity. This is 
precisely the situation in the decaying process in which the satellite 
makes thousands of orbits during its lifetime. The beauty of the 
technique is that it circumvents the problem of computing the 
trajectory for each orbit which can cause numerical integration 
schemes to diverge over the long integration intervals. 

In the variational equations (3-12) the radial distance will be 
written as 
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so that the exponential function 

exp r /3(r - r)l = exp[ /i{a -a- a e ) + fjaecosEj 

P ' o o o 

o (4-Z; 

where a and e are the initial semi-naaior axis and initial 
o o 

eccentricity respectively. To facilitate the theory a new dimension- 
less variable x will be used to replace e. 

x = |3ae (4-3) 


Differentiating Eq. (4-3) gives 
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Now it is evident that the x varies like e since x passes through 
stationary values when cos E = - e , but, on the average, decreases 
with time. 

The purpose of introducing the new variable x becomes clear 
th 

after considering the n order modified Bessel function, I^(x) . 
the first kind 

1 

I (x) = — / cos nE exp (x cos E) dE 

n V 

o 

along with the relation 

1 

0 = •— f sinnE exp(xcosE) dE . (4-5) 
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o 

The average equations from Eqs. (3-12) are 
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After setting w = 00 ^ , it is apparent that the average equations 
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can be written as power series in e, for small e, with Bessel 
Function coefficients as functions of the dimensionless variable x. 


4. 2 The Basic Equation for Orbit Contraction 


The most dramatic and important effect of the dissipative 
force of atmospheric drag on the satellite's orbital path is the 
reduction of both the semi-major axis and the eccentricity. In 
general, the orbit contracts from an initial elliptical orbit to a 
nearly circular orbit so that e and x approach zero in the final 
stages of the decay. At this time the satellite is in the last few 
orbits and entry is imminent. The main purpose of this chapter is 
to provide a rigorous mathematical treatment of the orbit contrac- 
tion problem. 

For small eccentricity , the integrand of the first equation of 
(4-6) can be expanded in a power series in e and integrated to 
obtain 
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Following the same procedure for the x equation 
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Let 


+ lla < 78 1^ + 3113+315) + 0(e^ 


(4-9) 


be the dimensionless semi-major axis. Dividing Eq. (4-7) by Eq. 
(4-8) and expanding the ratio in a power series 
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where the ratios of Bessel functions have been defined as 


n ^ 1 


(4-11) 


The Bessel functions satisfy the recurrence formula 




so that any y (x) can be expressed in terms of y (x) and x . 
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Putting the eccentricity, e, in terms of z and x 


e = e — 

z 


is a small quantity of order 10 or less. 

Thus, from Eqs. (4-10) - (4-15) the basic equation for orbit 


contraction is derived 
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The basic equation is a nonlinear differential equation with 
varying coefficients. The true nature of the equation was not 
recognized by King-Hele, who worked with Eq. (4-10) to the order 

3 ^ 

e . In the form of Eq. (4-16) it is appropriate to apply Poincare's 
method of small parameters. Since the parameter e is very small, 
to a higher order in e , the solution of the basic equation can be con- 
sidered to be the exact solution of Eq. (4-10) truncated to the order 
included. 
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4. 3 Integration by Poincare's Method of Small Parameters 

Poincare's method for integration (Poincare, I960) of a 
nonlinear differential equation containing a small parameter is a 
rigorous mathematical technique, proven to be convergent for small 
values of the parameter e . It has been used extensively in analytic 
work in celestial mechanics (Moulton, 1920). 

The method begins by assuming a solution of z in the form 
2 3 4 5 

z = z+ez+e z+e e z. + e z+... (4-17) 

o 1 /i i 4 5 

After substituting Eq. (4-17) into the basic equation , Eq. (4-l6), 

and equating coefficients of like powers in e , the following 

differential equations are obtained 
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with the initial conditions 

z (x ) = 1 , z (x ) = z (x ) = ••• = 0 . (4-19) 

O O 1 O 2 o 

The success of Poincare's method is not guaranteed; it depends on 
whether the integration of Eqs. (4-18) can be put in terms of known 
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functions. 


In analyzing Eqs. (4-18), an important rc^currcnco formula 
was discovered which greatly facilitates the integrations. 
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where n ^ 0 and p(x) is an arbitrary function. The formula is 
derived by using the well-known relation 
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and for n = 0 
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or 
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/p(x) y^' dx = / p(x) + if “ 


p(x) n 
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r p'(x) , n 

/ y dx . 

n o 


Rearranging the equation gives Eq. (4-20). 

Eqs. (4-18) can now be integrated with the initial conditions 
(4-19). The solution for is simply 
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(4-25) 


From Eq. (4-22), 
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Since z =1, the z_ equation becomes 
o 2 
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Integrating 
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Applying the recurrence formula (4-20) with p(x) = x and n = 1 , 
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so that z^ with the initial conditions is 
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(4-29) 
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The solution of z„ , z . and z^ are obtained in the same 
3 4 5 

manner, but the integrations are much more laborious. The z.(x) 
can be expressed in terms of two functions 
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The final solution is 
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The semi-major axis of the orbit under contraction is given 
as a function of the parameter x. 


— = 1 +e z (x) + e ^z (x) + e ^z (x) + e ^z (x) + e ^z (x) 

a 1 2 3 4 5 

° (4-32) 


The accuracy of this solution has been tested by numerically integra- 
ting the basic equation (4-1 6 ) and comparing with Eq. (4-32). The 
expansions used to obtain the basic nonlinear equation are only valid 
for small and moderate values of eccentricity. To be exact, expan- 
sions in elliptic motion apply to eccentricities which are less than 
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0. 663. Above this value, the series are no longer absolutely con- 
vergent. However, in the accuracy test it was found that the analytic 

solution was always greater than the numerical integration with a 

5 2 

maximum error of approximately e e / 5(1 - e ) for 0. 1 < e <0. 99. 

o o — o — 

It is interesting to note that even as e -*1, which is outside the 

o 

region of strict mathematical validity, the maximum error is less 

_3 

than 1/ 10 pr (of the order 10 ). The reason for this extended 

^o 

range is that as e -*■ 1 , a oo and the parameter e -»• 0 . This 
o o 

important result prompted investigation of the asymptotic expansion 
of the basic nonlinear equation and the development of a closed form 
theory for orbit contraction for large eccentricity which is presented 
in Section 4. 6. For small values of eccentricity the solution is 
extremely accurate. For example, when e^ = 0. 1 and 6 = 0. 008 , 
Eqs. (4-31) and (4-32) provide 7 digits of accuracy, while for the 
same case King-Hele's heuristic solution gives only 4 digits of 
accuracy. Thus, the present solution provides a major improve- 
ment in that it is much more accurate and that it is uniformly valid 
for all eccentricities, (0 < e < e ) . 
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4.4 The Parametric Solution of the Orbital Elements 


The solution for the semi-major axis is given as a function of 

X by Eqs. (4-32) and (4-31 ) and is plotted against x/ x^ in Fig . (4-1). 

Initially the solution is nearly linear with the slope in the figure 

approximately equal to e , but near the end as x and e approach 

o 

zero , it exhibits rapid decay. This explains the difficulty 
encountered by King-Hele and why he treated the problem in separate 
phases in his analytic integration. 

The eccentricity can be recovered from Fig. (4-1) using 
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(4-33) 


From the parametric solutions of the semi-major axis a and 
the eccentricity e it is easy to find expressions for the periapsis, the 
apoapsis and the period. 

The drop in scale height of the periapsis is obtained from 


r - r 
P P 

= sir -a(l-e)| =6a -Ba e + B ae-B a (1+ez +...) 

H Ip l o'^oo ‘o 1 

o 




or 


r - r 

Po P 2 3 4 

jj— = (x-x^)-(aj+ez2 + € + e + e 


(4-34) 


Similarly, the other quantities of interest can be written in 
parametric form. 
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The drop in scale height of the apoapsis is 
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The ratio of periapsis to initial periapsis is 
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The ratio of apoapsis to initial apoapsis is 


z + c X 
1 + e 


(4-37) 


The orbital period is 
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z (x) 


(4-38) 


Eqs, (4-34) - (4-38) are plotted in Figs. (4-2) - (4-6) as functions of 

x/x . The parameters used are the initial eccentricity e and the 
o o 

initial periapsis distance r , or equivalently the dimensionless 

^o 

small parameter l/6r . Then 

Po 
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which tends toward zero as e^ 1 . By using the three values 
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= 0.005, 0. 01 and 0.02, a wide range of periapsis heights 


are covered. In naaking these figures reference was made to the 

1 


dependency of the accuracy on the initial parameters e and 

° P 

as mentioned in the previous section. The figures are plotted 
only for those values which provide imperceptible deviation from the 
exact solution. 


4. 5 Explicit Formulas for the Orbital Elements 

The solutions for the orbital elements have been obtained in 
parametric form, as functions of x in the previous section. In this 
section the elements will be obtained explicitly in terms of the 
eccentricity. One approach might be to eliminate x between the 
parametric equations. However, because of the transcendental 
nature of the solutions, this task is cumbersome. Fortunately, 
since e is small, x can be eliminated by applying Lagrange's 
expansion. 

To derive the explicit expression for the semi-major axis 
in terms of the eccentricity, z must be written in a form suitable 
to Lagrange's expansion: 

z = p + e <|)(z) (4-40) 

and since 
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it can be written as 
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so that a is proportional to e. Then 
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Lagrange's expansion of Eq. (4-40) is 


^ , p + ^ ,p)] (4.45 

where the expansion is carried out first and then p is set equal to 1 
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to give to the order of e 
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In the expressions for h,(a) the following definitions were used 

a I (a ) 

A = ay (a) , z = Log . (4-48) 

o 1 X l(x ) 

o 1 o 

Thus, Eqs. (4-46)-(4-48) provide an explicit expression of the 

variation of a/ a as a function of e/ e since a - x (e/ e ) . 

o 000 

Following the approach of the previous section, the other 
orbital elements can be easily expressed in terms of e to obtain the 
drop in periapsis, the drop in apoapsis and the change in orbital 
period. 
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o 

The last expression, Eq. (4-53), provides a powerful method of veri- 
fying the assumption made on the atmosphere from two of the most 
accurately and easily measured orbital elements, namely, the period 
of revolution and the eccentricity. 

Eqs. (4-46) and (4-49) -(4-53) are plotted in Figs, (4-7)-(4-12) 
respectively. 

Fig. (4 - 7) plots the solution z = z{a) = z(e/e ) as a function 

of the eccentricity for several values of e^ , The range of validity is 

limited to 0 < e < 0. 5 since the z(a) solution is not as accurate as 
o — 

the z(x) solution. It was found that z{a) exceeds the numerical solu- 
tion (found by integrating the basic nonlinear equation (4-16)) by a 
maximum value approximated by ej^/15(l-e^), which still gives 7 
digits of accuracy for e = 0. 1 , but diverges as e -*■ 1 . For e = 

0. 5 the error is imperceptible in the figure. The decay in the 

periapsis distance r versus the eccentricity for e =0.1 and e = 

p o o 


44 



remains nearly 


0.4 is presented in Fig. (4-10). The ratio r / r 

equal to one for a large portion of the decay process, but as e -* 0 

the drop in periapsis altitude increases rapidly. For small values 
1 


of 




(0. 005), which correspond to higher initial periapsis 


altitudes for fixed p , decay is much slower than for large values 
(0. 02) as evident in the figure. The fractional error in this plot 
and in the others is kept to an imperceptible amount by considering 
the error formula mentioned earlier. 

The decay in the apoapsis distance r^ versus the eccentricity 
for several values of e^ is presented in Fig. (4-11). It is evident 
that the ratio r / r decreases rapidly with the eccentricity. 

SL 3, 

^ 1 

Initially the parameter seems to have little effect, but as the 

P ^p 
*^o 

eccentricity approaches zero the larger values of 


1 




yield 


more rapid decay. 

Finally, the decay in the orbital period T as a function of the 
eccentricity is presented in Fig. (4-12). As pointed out by King- 
Hele, this functional relationship T = T(e,e) provides a powerful 
formula for testing the atmospheric parameter e since the orbital 
period and eccentricity can be accurately measured. 
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4. 6 The Contraction of Highly Eccentric Orbits 


For the case of orbits with large eccentricities, King-Hele 
used an entirely different method from the following to derive the 
basic equation for orbit contraction. Using the present notation, his 
equation takes the following form 


dz 

dx 
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z + ex 


(4-54) 


This equation can be directly derived from the basic non- 
linear equation (4-1 6). Since x = p ae, when e -* 1 , a -► oo , x 
becomes very large and the asymptotic expansion for Bessel's ratio 
y (x) is 
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(4-55) 


Substituting Eqs. (4-55) into Eq. (4-l6) provides 
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z 2 " 3 4 


(4-56) 


which is immediately recognized as the development of Eq. (4-54). 
King-Hele provides an approximate solution to the nonlinear 
equation (4-54) by assuming that on the right-hand side z is 


approximated by z = 1 - e (x - x) and then integrating. In this 
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analysis, the exact solution will be given. 
Using the transformation 


z = e (x + q) (4-57) 

and changing the independent variable from x to q results in the 
Bernoulli equation 


dx _ 4x 
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dq q 


(4-58) 


Substituting the change of variable 
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into Eq. (4-58) provides 
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Upon integrating it is found that K can be expressed in terms of the 
exponential integral . 

2q 

K = - 4 r ^ d(2q) + C 
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= -4 E.(2q) + C . (4-61) 
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Thus, the exact solution in parametric form is 
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along with Eq. (4-57) and the initial conditions 


z(x ) = 1 
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% = 7 
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(4-63) 


The exponential integral can be evaluated through its 
asymptotic form since the argument 2q is large. In general, the 
exponential integral of the form 


E (x) 
n 

can be integrated by parts 
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By repeated application of this formula, the asymptotic expansion 
for large x is deduced 
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(4-65) 


When n = 0, Eq. (4-64) becomes the exponential integral appearing 
in Eq. (4-62). Taking 6 terms of the series (Eq. (4-65)), for 
X > 50, the solution is identical to numerical values tabulated by 
Abramowitz and Stegun (1972). 

Numerical computation has revealed that the analytic solution 
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z(x) of the basic nonlinear equation (4-l6) and the present solution 

z(q) of the asymptotic equation are in very close agreement. Even 

for e = 0. 99 the two solutions are nearly identical except for very 
o 

small values of x/ x^ . This gives further testimony to the wide 
application of the solution z(x) of the basic nonlinear equation. 


4. 7 The Contraction of Nearly Circular Orbits 


For very small eccentricities, x is very small and the series 

expansion of y (x) for small x is 
o 


2,1 13 1 5 

0 x 4 96 1536 


(4-66) 


The derivation is discussed in detail in Chapter 5. Substituting 
Eqs. (4-66) and (4-13) into Eq. (4-10) results in the following 


1 ^z 

e dx 




(4-67) 


Eq. (4-67) is recognized as the expansion of 


€ dx 


x(l+3- 

X 


(4-68) 


C X 

Writing e , the equation can be put in the form 

z 
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dx 


( 1 + 3 - ) dz = — 

Z X 

Upon integrating 

z + 3e Log z = 1 + 2e Log — (4-69) 

X 

o 

a closed form parametric solution is obtained for very small values 
of eccentricity. It would be useful to have the explicit solution of the 
semi-major axis z in terms of the eccentricity. This is accomplish- 
ed by using Lagrange's expansion as before, with 

0 

tj) (z) = 2 Log — - Log z 

Q 

O 

Applying the expansion, 

e 2 3 

z(e/ e ) = 1+6 2Log — (1-6+6 -6+-..) 

o 

= 1+7^ (4-70) 

o 

provides the closed form solution of z as a function of eccentricity 
for very small values of eccentricity. 

4. 8 The Contraction of Circular Orbits 

For contraction of circular orbits, it is necessary to go back 
to the variational equations (3-12) and substitute e = 0 into the 
equation. Since r = a , the equation becomes 
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dr 

dM 


2Z r 
o 


P(r -r) 


(4-71) 


where M is the mean anomaly. The equation can be rearranged to 


give 


p r 2Z P 

e d p r = - e ° dM 


(P 


rf 




(4-72) 


The first integral is simply E ^ (p r) , as defined by Eq. (4-64) 
with n = -1. Assuming initial conditions r and M the integration 

Po 

provides 

1 _ 

E , (p r) - E , (p r ) = - 2e Z e ^ [M - M 1 (4-73) 

-1 -I p o ‘ o 


The exact total time in circular orbit is found by setting the 
radial distance to the minimum final distance, r^ , at which the 
orbit can barely be maintained and by replacing the mean anomaly 
M by ^ ^ th e init ial mean anomaly by zero 




_f_ 

M- 


2e Z e 
o 


T7T 

n -J 


(4-74) 
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4. 9 The Orientation of the Orbit During Decay 


Due to the fact that oblateness has been ignored in the math 
model, it turns out that the argument of the periapsis, w , remains 
a constant 


w = w 


(4-75) 


from the third of Eqs. (4-6), This allows the replacement of sin a 

and cos a with terms in sin E and cos E in the last two variational 

equations for , the longitude of the ascending node, and i , the 

inclination of the orbit plane. This is done by using the well-known 

relations of true anomaly a - and eccentric anomaly E 

o 

cos E - e 

cos (a - caj ) = — 

o 1 - e cos E 


, , e sin E 

sin (a - OJ ) = ; = — 

o 1 - e cos E 


VT 


n 


(4-76) 


From Eqs. (4-76) it is easy to obtain 


cos O' = 


sin a = 


(cosoj )(cosE-e) - (sinco ) 
o o 




sin E 


1 - e cos E 


(sinco )(cosE-e) + (cosco ) \jl - e sin E 
o o » 


1 - e cos E 


(4-77) 


The terms sin a cos a and cos a are needed for the variational 
equations. In terms of eccentric anomaly, they are 
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sin a cos a = 


cos w s 
I o 


2 2 2 1 
imj (cosE-e) +(e -l)(l-cos E) I 


(1 -e cos E) 


+ (2 cos^w -l)j^-y^l-e^ sin E(cos E-e)j 


2 1 2 , ^ ,2 
cos a = !T“ <cos w (cos E-e) 


(1 - e cos E) 


'o -f-' 


- 2cosw^sin« \/l-e sinE(cosE-e) 


+ sin^w (1 -e^)(l -cos^E ) 
o 


(4-78) 


It is now clear that the average equations can be obtained in the 
same manner as in Section 4. 1 by using Eqs. (4-5). Noting that the 
sinE terms vanish, the average equations for SI and i become 
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dE 


- r wZ 

o 


5/2 


r===Tf=^ Vo'] 

v'Pf/'p V'-" Po 


1 

— — / COSO) sino) 

2ir o o 


X (1-ecosE)^^ (1+ecosE)^^^ |^(cos E-e)^+(e^-l)(l-cos^E)J 


X exp (xcosE) dE 
-r wZ 


1 di 


P o P 

= /==rF=?' ( r ) ="p p 

yu£/r Vi-« Po 


5/2 


sini dE r 


1/2 1/2 

X ~ /"(1-ecosE) (1+ecosE) 

2Tr 




r 2 , 2 2 2 2 1 

X cos oj (cos E-e) +sin co (1-e )(l-cos E) 

L o o J 


X exp (x cos E) dE 


(4-79) 


Expanding to the first order in eccentricity and identifying the 

Bessel functions, I (x), the average equations become; 

n 
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dE 


- r wZ 

Po ° 




a-a e 
o 



cosw sinco 
o o 


X [l^ - 2elj + 0(e^)] 


d(Log(tan -)) 
dE 


-r w Z 

P„ o 




5/2 


— I — - — 1 exp (3 (a -a-a e )| 

->1^/ [ o ooj 


X 


( 1 1 2 2 2 2 
I + — I^(cos Cl) -sin Cl) ) - 2ecos o) I +0(e ) 
l.2o22 o o ol ' 


(4-80) 


Dividing both equations of (4-80) by the average 


dx 

dE 


equation, 


Eq. (4-8) 


dP 

dx 


C w 



[l^-2el,+0(e‘')] 


T +1 e (31 +1, ) + 0(e^) 
1 Z o Z 
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d(Log(tan-)) 




2Vl-e 


[i 


1 2 2 2 2 
— I_(cos CO -sin CO ) - 2e cos co I, + 0(e )| 
2 2 o o o 1 


[ll +'><"'>] 


(4-81) 


Since the term e w is very small, and i vary slowly. For the 
present development, a first order analysis will be adequate. This 
implies setting e = 0 in Eqs. (4-81). 


dS2 e w ( o . . . 

"5 = -rr- J — r cos CO sinco y,(x) 

dx 2»pf o o 2 


d(Log(tan|)) aj 


V [i 1^0 + i o- o’l'z] 


Substituting y^ = y - “ into Eqs, (4-82) 
^ o 


dS2 e w / o . , 2 

— = -r — -t/ — — cos CO sinco (y - — ) 
dx 2 ^ pf o o o X 


d(Log(tan-)) 


2 2.21 
cos (o y -(cos CO -sin co ) ~ 
o o o ox 


(4-83) 
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Finally, Eqs. (4-83) can be integrated from to x with initial 


conditions and i. 


S2 = + 

o 


e w 




\ g.f 


cos (i) 


sin w z,(x)-2Log— 1 

o oL‘ 



gw 

2 



2 

cos CJ z 
o 1 


(x) 


2 2 

(cos o) - sin oj ) Log 
o o 



(4-84) 


where z^(x) is defined by Eq. (4-26), Thus, first order solutions 

have been found in parametric form for the longitude of the ascending 

node, , and the inclination angle, i. For small eccentricities, 

Eqs. (4-84) indicate that varies very slightly, not only because 

w is small, but because as x -► 0 , z, (x) -► 2 Log — . Similarly, 

1 X 

o 

for the inclination angle, it can be shown that for small x 


(» 6 1 

XI o , e w 

— 1 tan — , where 6 = 

X / 2 2 

o ' 

Since 6 is very small, the inclination of the orbit plane hardly 
changes during the lifetime of the satellite for nonzero values of x. 
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(Tp ”^p)/ ^ ~ drop of periapsis in scale heights, 


e = initial eccentricity and x = p ae = ae/H. 

T^^ie plot is generated from Eqs. (4-34) and (4-31). 
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(RflO“Rfl)/H 

70.00 105.00 140.00 175.00 



Fig. 4-3. Variations of (r_ - r_)/ H as Function of x/ . 

o 

(r -r^)/ H = drop of apoapsis in scale heights, 
e^ = initial eccentricity and x = pae = ae/ H. 

The plot is generated from Eqs. (4-35) and (4-31). 
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RP/RPO 



.o. 

e = initial eccentricity and x = pae . 

l^e plot is generated from Eqs. (4-36), (4-32) 

and (4- 31). 
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Rfi/RflO 



o 


- ratio of apoapsis to initial apoapsis, 

= initial eccentricity and x = pae. The plot 
is generated from Eqs. (4-37), (4-32) and (4-31). 
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x/xo 


Fig. 4-6. Variations of T/T^ as Function of x/xq . 

T/T^ = ratio of orbital period to initial orbital period, 
e^ = initial eccentricity and x = pae. The plot is 
generated from Eqs. (4-38), (4-32) and (4-31). 
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Fig. 4-7. Variations of z as Function of e . 


z = nondimensional semi-major axis = a/ a 
e = eccentricity. The explicit solution z(e) 
generated by Eqs. (4-46) and (4-47). 


and 

is 






Fig. 4-9. Variations of (r -r )/H as Function of e. 

do O' 

(r -Tg^)/ H = drop of apoapsis in scale heights 
“o 

and e = eccentricity. The plot is generated by 
Eqs. (4-50) and (4-47). 




as Function of e. 


Fig. 4- 


0. Variations of r / r 

P P. 


r^/ r = ratio of periapsis to initial periapsis 
^ ^o 

and e = eccentricity. The plot is generated by 


Eqs. (4-51), (4-46) and (4.47). 








T/TO 



Fig. 4-12. Variations of T/ as Function of e. • 

T/ = ratio of orbital period to initial orbital 
period and e = eccentricity. The plot is genera- 
ted by Eqs. (4-53), (4-46) and (4-47). 
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CHAPTER 5 


TIME IN ORBIT 


5. 1 The Time Equation 


The final orbital variable to be determined is the time which 


is fotind from Kepler's equation 


I 


a 


E - e sin E 


(5-1) 


In its differential form 


d t 
dE 


3/2 


The average equation is simply 

3/2 T 


( 1 - e cos E) 


(5-2) 


3/2 


dt 

dE 




- ( 4 ) 


(5-3) 


where T is the initial orbital period. Since the other elements 
o 

have been given in parametric form as functions of x, it would be 
most convenient to find t as a function of x also. Dividing Eq. (5-3) 
by the average equation for x, Eq. (4-8), gives the average equation 
for the time 
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dt 

dx 


(5-4) 


T r exp(x ) 
op o 

o 


4ir 6a z 
o o 


exp 


[pa^(z-l)] 


1/2 




The exponential term is found by using the parametric solution for 

z = a/ a from Eq. (4-32). 
o 


exp 


[f5ao(a-u] 


r 2 -| 

= exp z^+ez^+e z^ + . . . 


c 

= exp(Zj) l^l+e Z^+ -^(z^ +2z^) + , 


xl, (x) 


2 


rm [l+C», + ^(a,-+2z,) + 


X 1 (x ) 
o 1 o 


2 2 ' 2 




(5-5) 


Defining the dimensionless time t as 
2ir6a SfC^ 

° Pq ° 

T = — = X I (x ) exp (-X )t 

ml o 1 o o 


(5-6) 


the dimensionless time equation takes the form 


dr 

dx 




1/2 r , ex 

= L 2^ 


2 2 

( 3 y +7^) + ^-^(11+73) + 

8z 


• • • 


(5-7) 


After substituting for z(x) , applying the binomial expansion and 
dividing, Eq. (5-7) becomes 
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dr 

dx 


-X 1 - (z, -2z rf- 3xy +xy ) 

. ^ 1 u O u 


e 2 2 , 2 2 2 2 

+ -rr (-1 lx -X y_ + 18xy z, + 6xy_z, + 18x y + 2x y. 

8 3 ol 21 o 2 


2 2 
+ 12x y y_ - 12xy z - 4xy z + 3z - 4z 
o 2 o 2 2 2 X 

- 4»jZ^ + 4z2^ + 8:=3)] 


(5-8) 


5. 2 Integration of the Time Equation 


It is convenient, but not necessary, to assume the solution 
of the time equation in the following form 


T= T +eT +e T_ 
o 1 2 


(5-9) 


Substituting Eq. (5-9) into Eq. (5-8) provides the three differential 

equations for t , t, and t_ 
o 1 d 


dr 
o 

dx 

dr 

1 

dx 

!l2 

dx 



(2A -1) X + “ xz 
o 2 1 

-2x^ + -^(7x ^-8A ^-llA )x + 9x^y 
Zoo o o 

- (10 + 7A ) xz - xz ^ (5-10) 

O 1 o X 
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with initial conditions 


o o 1 o c. o 


(5-11) 


Integrating Eqs. (5-10) from to x , it is found that 

1.2 2 
T = - (x - X ) 

o 2 o 


^1 = +4^^^ 


1 2 
- 7 J ^ 


14 4 

= - ^ (X -X ) 

2 o 

1 2 2 2 2 
+ - (x -X )(7x -8A -llA ) 

4 o o o o 


12 63 2 2 

-xZj(10+7A^)-.j^x.j 


+ 1 [28 + 7A^] / kS^ dx 

X 


63 r c 
-— / z X y dx 

8 1 ^o 


(5-12) 


c 2 r 2 

Unfortunately, the two integrals J x y dx and J z x y dx cannot 


be expressed in terms of known functions. Approximate techniques 
will be applied in the next section to yield an accurate solution of 
these integrals. " 


73 



5. 3 Approximate Integration of the Unknown Integrals 


For the modified Bessel Function I (x) the following series 

n 

expansion is available for small x (Beyer, 1976) 


I (X) = -^ |l + 

" Z”n! ( 


2 4 

2 . 1! (n+1) 2 • 2! (n+l)(n+2) 


+ 


2 • 3 ! (n+l)(n+2)(n+3) 


4* • • * 


(5-13) 


For n = 0 and n = 1 this formula becomes 


T / \ 1 j. / ^ 1 yX,6, 1 ,x,8 . 1 ,x,10. 

I (x) = ! + (-) +Tt(t) +T7T(t) + 


2' 4 '2' 36 '2' 576 '2' 14400 '2 


(t) + 


T / \ - 1 ,X,5. 1 ,x,7. 1 x,9 . 1 , X 11 

^ 1 ^^^ 2 '^ 2 ^ 2 ^ "^ 12 ^ 2 ^ ■^ 144 ^ 2 ^ "^2880 ^ 2 ^ '^ 86400 ^ 2 ^ 


(5-14) 


Dividing the first equation of Eqs. (5-14) by the second gives the 
series expansion of y^(x) for small x which will be designated y^ 


2.1 13.1 5 1 7, 13 9 

^o " X 4^ " 96^ 1536 ^ " 23040 ^ "^4423680 ^ 

s 


(5-15) 


On the other hand, for large x, the asymptotic expansion 
(Abramowitz and Stegun, 1972) is 
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I (x) = 
n 


■^'ZlT X * 


t^-1 + (t^-U(tx-9) 


8x 


2! (8x) 


(tx-l)(ti-9)(tx-25) 

3!(8x)^ 


(5-16) 


where |jl = 4n . Setting n = 0 and n = 1 gives the asymptotic 


expansions for I and I. . 

o 1 


I (x) = - 

° V 


2ir X 


- + ~ 
T I 8j 


9 . 75 . 3675 

+ — + — + ^ 

^ 128x 1024x 32768X 


59535 


262144X* 


+ • • • 


X 


I,(x) = 


V 


2ir X 


15 


105 


4725 


128x^ 1024x^ 


32768X 


72765 

262144x' 


(5-17) 


After performing the division , the expansion of y^(^) large x 


designated by y is 
°L 

1 3 3 . 63 , 27 

y = 1 + — + j + 7 + 7 + 

°L 8x 8x 128x^ 


32x~ 


+ . . 


(5-18) 


Multiplying Eqs. (5-15) and (5-18) by x and integrating term by 


term from x to x , the approximate solution of the first unknown 
o 
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integral is obtained for small and large values of x. 


X 


2 ^ ,2 2^ 1,4 4, 1,6 6^ 

X y dx = (x -X ) + — (x -X ) - •' / (x -x ) 

^o o ' 16 ' o' 576 ' o ' 

X s 

o 


/ 


+ ^— (x®-x ®) - — i— 

12288 ' o ' 230400' o ' 


13 , 12 12, 

53084160 ’^o ^ 


for X < 3 


"" 2 

X y dx 

X °L 

o 


/ 


1,3 3, 1 , 2 2, 3 

— (x -X ) + -(X -X ) + - (x-x ) 
3 o 4 o 8 o 


^ 3 , X ^ 63 ' 1 1 . 27 . 1 I 

O o XX 

o 


for X > 3 


(5-19) 


Comparison of the approximate solution (5-19) and the numerical 
2 

integration of x y indicates a maximum error of 0. 3% . Since 
the integral appears in the e term, this is equivalent to an error 
of 0.003% or less, which is sufficiently accurate. 

Next, the second unknown integral which appears in the 
equation of Eqs. (5-12) must be evaluated. 

Integration by parts reveals 
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/ z^x^y^dx = /x^^dx- / |/x^y^dx jy^dx (5-20) 

X X L J 


The first integral on the right hand side has been solved and the 
second can be determined easily from the earlier work. By dropping 
the constants from Eq. (5-19) and then multiplying by Eqs. (5-15) 


and (5-18) and integrating term by term from x^ to x , the integration 
of Eq. (5-20) can be performed for small and large values of x. 


C* 2 , 1 4 1 

X s l_ 


1 6 1 8 1 10 13 12 

576 ^ 12288 ^ "230400^ '^53084160^ 


- (x^-x ^) - — (x"^-x "^) - — ^ — (x^-x ^) 
' o ’ 32' o ’ 3456' o ^ 


5 8 8. 299 10 10^ 

^ 147456 "^o ^ ■ 110592000^"^ “^o ^ 


623 . 12 12) 

+ 3r85049600 <^ '^o for x < 3 


r z X y dx = z 
‘'1 o^ 

X L 


[ 13 12 3 ,3 ^ 

JX +JX +-x+-Logx- 


63 27 

128x ~ 2 

64x - 


1.4 4. 5 , 3 3. 5 , 2 2. 

T2<* '-36<* -*o >-T6<* -*o ' 


1 . 9 , 1 1 ^ 9,1 1 , 

32 o 512 X x 256 2 2 

O XX 

o 

T r 3 3 3 , 9 9 1 

+ Logx --x+^ -32^08* + ^ + -2 
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- ’‘o L‘ t 

>— o 128x -I 

o 

for X > 3 (5-21) 

2 

The exact integration of z^x was determined numerically by 

computer and compared to Eqs. (5-21). The greatest error was 

2 

found to be less than 1%. Since the integral appears in the e term 
the error is of the order 0, 0001%. Thus, Eqs. (5-12), (5-19) and 
(5-21) give the solution of the dimensionless time t as a function of 
X to a high degree of accuracy. Now, for any value of the parameter 
X, the semi-major axis, eccentricity, periapsis, apoapsis, period 
and the time can be computed. The solutions are very useful in 
determining the total lifetime of the satellite as well. For example, 
the final value of one of the orbital elements such as periapsis or 
period can be specified for orbit decay and the corresponding x^ 
can be put into the time solution to obtain an accurate estimate of the 
satellite's lifetime. 

In order to provide a formula for the maximum lifetime, the 

limit of T , T, and T_ will be taken as x, -► 0 . This can be done 
o ’ 1 2 f 

by inspection for most of the terms in Eqs. (5-12). For the 
remaining terms, L' Hospital's Rule is applied. For example. 
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Llim 

x-*0 


Llim 

x-*0 


2 


X z. 



Lim 

x-*0 -2x 


= Lim 
x -^0 


1 3 

T X 



J_ 

96 


3 

X + . . . 


= 0 


(5-22) 


Similarly, 


2 2 

Lim X z = Lim xz Lim xz 
x-»0 x-»0 x-^0 ^ 


and 


Lim xz 
x-*0 


1 


Lim ■— 

x -^0 - 

X 


= Lim 


x-*-0 -X 


-2 


2 2 1 

= Lim - X ( - + T X 
x^O ^ ^ 


... ) 


= 0 


so that 

2 2 

Lim X z =0 (5-23) 

x-0 ^ 


Thus , the maximum lifetime of the satellite is the finite 
quantity 

max 
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2f 1 4 1 2 2 2 

+ e 4 X -t(7x -8A -llA ) X 

|_2 o 4 o o o o 

+ (28+ 7 A ) f x^y dx 

2 o o 

X 

o 


X 

o 


where the integrals are evaluated by Eqs. (5-19) and (5-21). It 
should be noted that the lifetime obtained from Eq. (5-24) is the 
maximum and that it is possible that the entry phase may commence 
at an earlier time, when x^ is not yet equal to zero, as mentioned 
earlier. 

The accuracy of Eqs. (5-24) and (5-12) along with Eqs. (5-19) 

and (5-21) has been tested by comparing the analytic solution to the 

numerical integration of Eq. (5-8). The greatest error found is in 

the fourth or fifth digit which is sufficiently accurate when the range 

of validity is restricted to small eccentricities, 0 < e < 0. 2. 

Since the maximum lifetime is expressed as a function of 

X and e , T, = (x ,e), in Eq. (5-24) the final equation can 
o Li E o 

max 

also be put in terms of e and e , = T- (e , e )• It is 

O Li Li O 

max 
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interesting to note that for very small values of eccentricity, the 
classical parabolic law is deduced. 


e 


2 


max 



(5-25) 


So that, as a rough estimate, the time in orbit is proportional to 

the square of the eccentricity. The expression for given 

max 

by Eqs. (5-24), (5-19) and (5-21) provides a substantial improve- 
ment both in accuracy and range of application over the classical 
parabolic law. 
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CHAPTER 6 


BALLISTIC ENTRY 


6. 1 The Ballistic Entry Equations 

Near the final stages of orbit decay the orbit circularizes 
and the entry phase is imminent. The universal entry equations 
(2-18) still apply. For C. = 0 and pr = constant, the equations 

J_l 


become 


dZ 

ds 

du 

ds 


= - p rZ tan y 


fs = 


2 >/^ 
cos y 

2 

cos y 
u 


Zu - utany 


dg 

ds 


= 1 


dn 

ds 


= 0 


f- = 0 


( 6 - 1 ) 


From the last three equations of (6-1) it is apparent that s 
is equivalent to a , the angle between the ascending node and the 
position vector, and that and i are constants. This reduces the 
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( 6 - 2 ) 


(6-3) 


dZ 

dc^ 

= - p r Z tan y 




2"^P r' Zv 

i’ y.. 1 

z 

fl 

dof 

cos y 

1 ■ v) 

dff 

= 1 - i 

V 


(6-4) 


Eqs. (6-4) apply to the motion of any ballistic vehicle in a 
nonrotating atmosphere from the initial time in orbit to the instant 
of contact with the surface of the planet. The atmospheric density 
varies exponentially with altitude and the p r term is assumed to be 
a constant approximately equal to 900. The dependent variables 
Z, V and y are the modified Chapman variable for altitude, the 
non dimensional velocity squared and the flight path angle. The 
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significance of the modified Chapman variable Z is that it allows 
a single trajectory solution for a given initial velocity and flight 
path angle to be computed which applies to every possible ballistic 
satellite, regardless of the mass, surface area or drag coefficient. 
This means that for particular initial conditions on the velocity and 
flight path angle, every satellite follows the same trajectory in terms 
of Z. The only difference is the particular altitude, which depends 
directly on the satellite parameters. Thus a single analytic solu- 
tion has a wide range of application. Furthermore, using the 
modified Chapman variable the resulting equations are free of the 
restrictions of the original Chapman formulation, as discussed by 
Vinh and Brace, 1974, so that the exact tables of Z functions for 
entry analyses can now be generated to replace those of Chapman 
and Kapphahn, 1961. 

In the next three sections, three distinct analytic theories 
are developed to cover the cases of ballistic entry from circular 
orbit, ballistic entry from nearly circular orbit and ballistic entry 
at moderate and large initial flight path angles. This is due to the 
fact that singularities arise in the treatment of the problem as a 
result of dividing either the first or the second of Eqs. (6-4) by the 
other two and integrating. Unfortunately it does not seem possible 
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to have a single analytic solution which is uniformly valid for all 
values of initial flight path angles because of the nature of the 
problem. In the case of atmospheric entry from circular orbit, the 
magnitude of the flight path angle, initially zero, rapidly increases, 
approaching 90 as the velocity becomes small. On the other hand, 
for steep angle entry, the flight path angle changes very little- -of 
the order of tenths of a degree--as the nondimensional velocity 
decreases from unity to one tenth the original value (between Mach 
2 and 3). The zero order terms of the two theories are totally 
different, which indicates that separate theories are necessary. The 
entry theory for large and moderate flight path angles has as its 
small parameter 

- 1 

^ ~ - 2 

3 r V. tan y. 

1 

where v. and y. are the initial nondimensional velocity squared and the 
initial flight path angle, respectively. Obviously, when the magnitude 
of y. becomes small a new small parameter must be found. For 
small initial flight path angles a simple power series solution is 
found by assuming the magnitude of is small and is zero . The 
same approach is not available for the zero initial flight path angle 
analysis, since both initial derivatives are zero and all the terms of 
the power series would be identically zero. Following the choice of 
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variables used by Yaroshevskii (1964), the problem is resolved in 
such a way that the initial point is a singularity, but has the proper 
behavior as the independent variable approaches zero. The 
Yaroshevskii solution is found as the zero order term and then a 
first order term is added, which greatly increases the accuracy. 

In the next sections the entry from circular orbit is 
developed first, since it represents the final stages of the majority 
of decaying orbits. Next, the small angle theory is presented and 
finally the moderate and large angle theory, which applies to the 

case in which entry is usually accomplished on the first orbit from 

( 


an initially highly eccentric orbit. 


6. 2 Entry from Circular Orbit 

The planar entry equations (6-4) will now be considered for 
entry from circular orbit. Dividing the first and third equation of 
Eqs. (6-4) by the second gives 
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dZ 

dv 


dv 

with Z and y as the dependent variables which are initially nearly 
zero and v as the independent variable, initially one and assunaed 
to decrease to less than one hundredth as the final condition. 
Precise initial and final conditions cannot be put on Z without the 
knowledge of the satellite parameters. To surmount this difficulty, 
the final orbit was defined as that orbit in which the initial Z was 
such that the satellite's velocity reduced to one tenth the circular 
velocity, while the longitude reached 360 degrees. In most cases, 
this definition provides the entire trajectory including or beyond 
the point of surface contact. Only in the case of very low density 
entry vehicles, such as balloons, would the definition have to be 
modified. 

The development of the theory is facilitated by the following 
change of variables: 
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X 


Y 

$ 

X 

e 


= 2Z 

= -yJF^ sinv 

= - Log V 
1 

pr 


( 6 - 6 ) 


Note that all the new variables are positive. Rewriting Eqs. (6-5) 
in terms of these, 


dY 

dX 


d$ 

dX 


$ 


1+^(1- Ze^) 


(l-e$^)(e^-l) 


Y 1^1 + ^(1 -2e^)] 
Expanding Eqs. (6-7) to one term in e 


dY 

dX 


$ 


1^1 + ( 2 e^-l)J 


2 X 

d^ _ (1-g ^ )(e - 
dX ■ Y 


— [l + ^ (2e^-l)J 


(6-7) 


( 6 - 8 ) 


Poincare's method of small parameters for integration of a system of 
equations can now be demonstrated. Assume to the first order a 
solution of the form 


$ 

Y 


$ 

o 

Y 


o 


+ e 


+ e Y, 


(6-9) 
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Substituting Eqs. (6-9) into Eqs. (6-8) results in two systems of two 
first order differential equations. 


dY 

o 

dX 

d$ 

o 

dX 


= $ 


X 

e - 1 


( 6 - 10 ) 


dX 


(2e^-l) 


d#. 


dX 


(e^-1) 


$ Y 

(2e^-l) - $ 2 1 

i o Y 

o o 


(6-11) 


Eqs. (6-10) must be solved first. They can be put in the fo 
of a second order differential equation. 


rm 


d^Y 


dX 


X 

e - 1 


( 6 - 12 ) 


with initial conditions 


Y^(0) = 0 , $ (0) = 0 

o o 


(6-13) 


An attempt to solve Eq. (6-12) by assuming a solution of the form 

2 

Y = a + aX + a_X + . . . fails because the initial conditions 
o o 1 2 

(6-13) force the a^^ to be identically zero. Apparently, the first 
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term is neither a constant nor a linear function of X. Yaroshevskii's 
approach will be employed to find the first term of the series. 

For very small values of X, v -► 1 and Eq. (6-12) becomes 




■° dX^ 


X 


(6-14) 


The solution of Eq. (6-14) is easily found by assuming Y = c X 


n 


and solving for c and n to obtain 


i. ^3/2 




dY 


dX 




(6-15) 


which satisfy the initial conditions (6-13) as X -»• 0 . 

With the insight gained from the analysis of Eq. (6-14), the 
solution to Eq. (6-12) will be 


with 


Y = y 


2 3 4 

y — 1 + a X + a X + a X + a X + . 

^ 12 3 4 


(6-16) 


(6-17) 


Using Eqs. (6-16) and (6-17) in Eq. (6-12) provides 


Xy + 4X 


2 


dX 





(6-18) 
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Expanding the exponential term and substituting for y and its deriva- 
tives allows the a. terms to be found after equating coefficients of 
like powers in X . Thus, 


Y = i 
° -v/T 


IX, 1 , X : 

^1 + 3 ( “ ) + 5 ( - ) 


47 

594 ' 4 ' 605880 


20021 ^ 4 

' 4 ’ 


'] 


(6-19) 


$ 


rp ^l/Z r, 5,X, 7 X 2 47 , X,; 


20021 . 2L \ ^ 

165240 ^ 4 ’ 


] 


( 6 - 20 ) 


Next, the second system of differential equations, Eqs. (6-11) 
corresponding to the € term will be solved. The resulting second 
order equation in is 


2 

d Y, 


dX 




(, X ,, r 3(2e^-l) 

— (e -1) Y ® 

O V o 


o] 


$ 

+ 2e^$ - (2e^-l) 

o Y 


( 6 - 21 ) 


The correct form for the series solution of Y^ can be obtained by 


n 


letting Y^ = c X and solving Eq. (6-21) for the lowest power of 


n 


X . From this it is found that the Y^ series solution must be in 
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the form 


^^J3 3/2 

■■ *' u ’ ' ‘ ^ V 

3 


= 1 + X + + 


( 6 - 22 ) 


Before substituting Eqs. (6-22) into Eq. (6-21), it is convenient 


to let 


C 

D 



Y 

o 


$ 

o 




(6-23) 


Now, with Eqs. (6-22) and (6-23) put into Eq. (6-21) 


2/3 M „ 2 


(4>'l 




+ 3 X 


+ X 


dX 


dX 


M X ^ r' 

2 ) + 4 "^>^1 = 


14 X ( 2 


D(C-B^) + Y B A (D + 3) 


j 


(6-24) 


Expanding A, B, C and D, 


, , 1 , X , , 1 , X ,2 47 , X ,3 , 20021 , X , 4 

^■^3^4^'^6^4^'''594 ^4^ 605880 ^ 4 ^ 


„ , ,5 , X, , 7 ,X ,2 47 , X ,3 20021 , X ,4 

^ 9 4 ^ 18 ^ 4 ^ 198 ^ 4 ^ ‘^165240 ^ 4 ^ 
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^ ,X, 8 , X ,2 8 , X,3 32 , X,4 

C = 1+2 (-)+ 3(-) +-j(-) + yg(^) 


^ , X, , , ,X,2 64 , X, 3 64 , X ,4 ,, 

D = 1 + 8 ( -) + 16 (~) + — { -) + — ( _ ) (6.25) 

After substituting Eqs. (6-25) and (6-22) into Eq. (6-24) and 
equating coefficients of like powers of X, the solution is found. 




tVT „3/2 
■■ ■ .X. 


, . A5 X ^ . 105047 X ^2 
^ 63 ^ 4 ^ 79002 ^ 4 ^ 


^ 191876677 ^21)^ 


132960366 ' 4 


(6-26) 


The solution for is 


^1 = 




I 189 4 


10^ X2 

f oor»rr> ' a • 


33858 ' 4 


191876677 
44320122 ^ 




(2e^-l) 


$ 


(6-27) 


Thus, the solution for reentry from circular orbit has been 
found in Eqs, (6-9), (6-19), (6-20), (6-26) and (6-27). It is 
interesting to note that the Yaroshevskii solution appears as the 
zero order term. With the first order term included, the accuracy 
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is greatly improved. The only drawback is that the expression, 

2 

Eq. (6-2 7), has a singularity at X = 0 due to the $ / Y term. 

o o 

However, the Yaroshevskii solution is not accurate near values of 
X cJi 0, and so the defect in the present theory is more than 
compensated for by the large gain in accuracy. 

6. 3 Entry from Nearly Circular Orbits 

If the entry phase initiates before the orbit has completed 
circularization, then a new theory must be applied which takes into 
account a small initial flight path angle. An approximate theory 
can easily be developed assuming 

sin Y (6-28) 

where the term on the left hand side is the deceleration in g's due 
to the atmospheric drag. This assumption is valid for a large 
portion of the entry trajectory. It is violated initially when the 
deceleration is very small, but this effect is only transient. The 
only other violation occurs at the final stage of entry when the 
magnitude of the flight path angle, y > becomes large and the 
dimensionless velocity squared, v , is small. By this time the 
essential features of the entry dynamics such as peak deceleration 
and peak heating rate have already occurred. With these consider- 


i 1 2 - V 

r Z V » ( — ^ — ) 
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ations, assumption (6-28) is justified and Eqs. (6-4) become 


dZ 

dof 


- P r Z tan \ 


dv 

da 


2 ^p r' Z V 
cos "Y 


lx 

da 


1 - 


V 


(6-29) 


Dividing the first and third equations of Eqs. (6-29) by the second 
equation , 


dZ 
d V 


2v 


siny 


lx - (1 - v) cos y 



(6-30) 


Now, Eqs. (6-30) can be transformed by the change of variables 


Y 


2 Z 


$ = - siny 


X = - Log V 

e = — (6-31) 

pr 


to 
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dY = ^ 

dX 

di ^ e ( e^- 1)( 1 - i 

dX Y 

with the initial conditions 


(6-32) 


Y(0) = 0 

$(0) = “ - sin-y. (6-33) 

The equations will be integrated by Poincare's method of 
small parameters. Assume 

Y = Y + e Y, + e ^Y_ 

o 12 

$ = $ +e$, + e^$- (6-34) 

o 1 2 ' ' 


Substituting Eqs. (6-34) into Eqs. (6-32) gives the three systems 
of two first order differential equations. 


dY 
o 

dX 


$ 

o 


d $ 

2. = 0 (6-35) 

dX 
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dY, 


dX 




d 


dX 




d Y, 


dX 



(6-3(.) 


dX 



- 2 
$ 

o 



- 2 $ $ 
o 


(6-37) 


The initial conditions are 


$ (0) = 

$ 

- sin V. 


o 

o 

1 


Y (0) = 

o 

Yj(0) 

= ^^(0) = 0 


ij(0) = 

i,(0) 

= 0 

(6-38) 


The solution of Eqs. (6-35) is simply 

i (X) = i 

o o 

Y (X) = i X (6-39) 

o o 


To solve Eqs. (6-36), first expand the exponential term and 
find . 
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So that 



(6-40) 


(6-41) 


Next, 


can be found by integrating Eq. (6-41) . 



(6-42) 


In a similar way the last set of equations (6-37) can be integrated to 


give 


^^(X) = - 


(^)( 5 . 


- 2 2 
(1 - $ ) 
o 


- 3 
$ 

o 


oo 


E 

m=l 


S 

n=l 


(n+m)m! n(n+l) ! 
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( 1 - $ % „m+n+l 

o \ X 

— I m=l n=l (m+n+l)m*m! n*n! 


, — 2 

00 

CO 

1 1 - $ 



0 . . 


E - 

— 

1 ri^l 

^1 ( 

1 ^ 

0 


- 2 2 
(l-$ ) 

00 

00 

0 

E 

E 

— 3 

$ 

^0 

m=l 

n=l 


(l-$ ) “ “ „n+m+l 

o v' X 

m=l _i (n+m+l)(n+m)m! n(n+l) ! 

®o ^ 

(6-43) 

The particular form of the solutions makes computation by computer 
very easy. It would be convenient to have the first few terms of the 


Explicitly, 






~ ) /ix2+ix^ + -l-x^+-!-x5+ 

- 3 I 4 ’' ^9^ ^144^ + 800 ^ +• 

m » 


144 800 


Y^CX) = 




o 

- 2 \ 2 


(ilfs-) Ux%^ 

-3 (l2^ +36 


x\-y-x^+-^ 

144 4800 


h • •) 


(6-44) 

Thus, the entry from nearly circular orbit problem has 
been solved to second order in Eqs. (6-34), (6-39), (6-41) , (6-42) 
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and (6-43). The appearance of the(l-<& $ term indicates the 

o o 

limitation of the theory to small $ ; very small values and large 

values of $ are excluded. For very small $ , the theory of 

o o 

entry from circular orbit can be used. For large and moderate 
values, a third analytic theory must be developed. With these 
considerations in mind, the theory for small initial flight path angles 
is restricted to the range 2. 5° " "V- • 

Next, the solution will be used to determine the speed at 
which the peak deceleration occurs. The deceleration in g'sis 
given by 

G = -| p r Y(X) e‘^ (6-45) 


Differentiating Eq. (6-45) with respect to X and setting the result 
equal to zero gives the equation for X at which maximum decelera- 
tion occurs. 


dY(X) 

dX 


Y(X) 


(6-46) 


Putting the expressions for the Y. into Eq. (6-46) gives the formula 
for the critical value of X and consequently the speed during 
maximum deceleration. 
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oo 


X 


1 + 


2 

tan Y 


(n+l-X) x” 
n* (n+1) ! 

n=l 


2 

£ 

■ 


2 oo 

i T'y 

CO 

^m+n+1 

2 

tan Y. 
1 



1 

n=l 

(m+n+l)m- m! n- n! 


4 

tan Yj 


oo oo 


,n+m 


(n+m+1 - X) X 

Zy Z/ (n+m+l)(n+m)m! n- (n+1) ! 
m= 1 n= 1 


(6-47) 


As a crude approximation, X = L for maximum deceleration which 
gives V = e ^ = 0. 367 88 or 


V 




- 1/2 

= e =0, 60653 


(6-48) 


This is the classical solution for steep angle entry in which the flight 
path angle remains nearly constant and the peak drag force occurs 
when the speed has reduced to 0. 60653 of the circular speed. Using 
X = 1 in Eq. (6-45) gives a rough estimate of the peak g 's. 


G 

max 



-1 

r siny^ e 


(6-49) 


Inspection of Eq. (6-47) reveals that for small angle entry 
the critical value of X is larger than 1, which corresponds to lower 
speeds. One method of obtaining the correct value of X would be 
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by numerical computation. Since X = f(X) , the first value of X must 
be guessed (X = 1 would be a good start) and a second value is 
generated from the equation. If jf'(X)| < 1 , which is generally 
valid here, the sequence X. = f(X ) , X_ = f(XJ, ... , converges 

X O w J. 

to the unique root X* = f(X* ) . 

It would be most convenient to have an accurate explicit 
solution for the critical value of X. Fortiinately, Lagrange's expan- 
sion can be applied (Eq. (4-45)) to Eq. (6-47) up to the first order 
term to give a good approximate solution. 


X 



tan y. 
1 



(n+l-X) x” 
n(n+l) ! 


(6-50) 



This simple formula gives the explicit value of X and, hence the 
speed, at which the peak deceleration occurs. For small initial 
flight path angles, it correctly indicates that the peak drag force 
occurs at lower velocities. In the range of small initial entry angles. 
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2,5° < - y. < 7° , the expression accurately predicts the critical 

speed to within 2% of the actual value. 

6. 4 Ballistic Entry at Moderate and Large Initial Flight Path Angles 

Although it is unlikely that a satellite undergoing orbit decay 
due to atmospheric drag would begin entry phase at a large flight path 
angle, there are certain cases such as return from the moon and 
atmospheric capture of probes to other planets in which an analytic 
theory for large initial flight path angles would be highly desirable 
and useful. In order to give a complete description of the entry phase 
for ballistic trajectories, the large and moderate initial flight path 
angle theory will be developed. 

Eqs. (6-4) will once again be used, but this time the first 
equation will be divided into the other two, providing 


dv 

fl 

siny . , 2 1 


dZ 

siny L 

2^^ Z ^ J 


.lx = 

dZ 

1 - V 


(6-52) 


p r Z vtany 




Eqs. (6-52) are ideal for analyzing the entry phase at steep entry 
angles. The dependent variables are the nondimensional velocity 
squared and the flight path angle y , while the independent variable 
is now the modified Chapman variable Z. When v^ = 1 the flight path 
angle remains nearly constant throughout most of the entry phase 
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while Z rapidly increases and v decreases. Initially, 


0 and 


lx = 

dZ 

remains small until the very end of the trajectory when v is small 
and the flight path angle suddenly decreases to nearly -90°. The 
known behavior of Eqs. (6-52) suggests an appropriate change of 
variables. 

Let the variables be transformed according to the following 

sin y. 

S" = 

sxn Y 


so that 


2Z 


q - - 


sin y. 


€ = ~ — = constant 

P r 


(6-53) 


d V 
dr| 


+ 


vS ( 






1^ ^ e ( V- 1 ) S 

dr| vp 



sin y. 


. 2 

sin y. 
1 


(6-54) 


Examination of Eqs. (6-54) reveals that Poincare 
parameters can once again be employed. 
Assuming 

2 


v=v +ev, + ev_ + ... 

o 1 2 

S = S +eS+e S + . . . 

o 1 2 


s method of small 


(6-55) 
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with initial conditions 


Since 


Eqs. 

three 


V = V 
o 


<V 


V. 

1 


S(ti.) = S. = 1 


v,(t|.)= V (q.) = 

1 1 Z 1 


= 0 


dS 
c 

dTi 


S (ti.) = S (t|.) = ... = 0 

1 1 2 1 


0 , S = 1 and Eqs. (6-55) become 
o 


(6-56) 


V = V +£V+€- v^ + .. 
o 1 2 


S- l+£S +£ S + . . . 

1 w 


(6-57) 


(6-57) will now be substituted into Eqs. (6-54) to obtain the 
systems 


dv 


dT| 


+ V =0 
o 


dS 
o 

dll 


(6-58) 


dv. 


dq 


+ V, + V S + 
1 o 1 


o 

•n 


2 

•n 


dT| 


2 

tan y. 


ri V r| ' 

. o 


(6-59) 


105 



0 


— + V + V S + V s + — = 

dT| Z 11 O Z T| 


dS. 


d T) 


— ( 1 - — + — j— 

. Z V T| . z 

sin y. o tan y. 


S S 
1 1 
n “v^n 


T) V 


(6-60) 


Note that the solution for has already been substituted into Eqs. 
(6-59) and (6-60). The integration of Eqs. (6-58) is easy 


V = 


S = 


-T| 

V. e 
1 


( 6 - 61 ) 


where 


V. 

1 


V. e 
1 


9. 

1 


(6-62) 


The zero order term is the classical approximate solution of the 
Chapman entry equation (Chapman, 1959). Note that any initial 
velocity and altitude are incorporated into the solution. 

Substituting the v expression into the flight path angle 
equation of Eqs. (6-59) provides 


dS 


1 


dp 


-A.- (1 --i 

tan\. ^ 


e^) 


(6-63) 


Integrating 
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^ = 


tan Y 


r- - - ■■■ 2— [Eih)-E;(v] 

i V. 1 V. tan V. *“ -* 

1 11 


(6-64) 


where £^(11) is the exponential integral 


r e^ 


(6-65) 


which is discussed in Section 4. 6. 


Now the equation of Eqs. (6-59) can be solved by variation 


of parameters. 


= C^(q) e 


- h 


( 6 - 66 ) 


Putting Eq. (6-66) into the v^ equation gives the equation for 


dC, 


= e 


B ^ ■ f] 


■y r\ V. 

2e 1 ^ 

- r Log 

■P t ^ h 

tan V. ' 

1 


IL 

i 


+ 


2 

tan y 


[EiC') - 


V. 


(6-67) 


Upon integrating 
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= 


2 

tan Y 


^ +'H) 


- V. Log ( n + tan Y-) 

I T). 1 

1 


— 'Hi 

+ v^(ti - Ti J - ( e^ - e ) 


( 6 - 68 ) 


The integrations of Eqs. (6-60) are performed in much the 
same way, but they are much more laborious. It is convenient to 
put the final results in the form 


V = V. 




(p) + e + 


...] 


s = 1 + c g^(p) + € ^g2(n) + ... 


( 6 - 69 ) 


where 


and 


e 


P 


S 


V. 

1 


1 


— 2 
0 r V. tan y. 
^ 1 1 



sin Y. 

^ 

sin Y 

P. 

1 

= V. e 
1 


(6-70) 
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V. Li - E 
1 o 

— ri — T) — 2 3 2 2 

v.(e -e ) -V. e L + v. (■r + tan v.) L 

1 1 1 2 1 


r — . 

I V. - V. T|. + e 

L 1 IX 


Tl 2 

+ e - 3 (1 + tan y 


.) l"| E 

X J o 


3 2 2 2 

+ (— + 2 tan V.) E - 2E (2t|) + v. tan v. F 

2 1 o 0 X 1 

n. , , 

— ~ tl X — - ^ ^ 

v.(ti - T|^) - ( e - e ) - v^(tan y- + 'H) (2 tan y^^ +"n ) 

[ Tl. 2 ~1 1 2 2 

e ^-v:(3 + 2tan yM ( t] - t|^) + -v^ (ri-ri.) 

— /3 w \ 1 

+ v^(3 -Tl )(e - e ) - - (e - e ) 


+ V 


- r— 2 2 

i 1 ViTii(3 + 2 tan y^) - (2 + tan y^^) e ^ 


— 2 Tl " 

+ v.(3 + tan Yj “'H )("n -Tl J + ("n - 2) (e - e ) 


1 — 2 4 2 2 

+ — V. (tan Y--^ti+ti ) L 
M X X 




2 _ _ 2 — 

+ I v.(2 + tan Yj) - 4 v^ti + v.ri + v.ti{3 -ti ) L + (2-r| )e 




+ 'tT1(ti-3)E ^ + 2(ti- 1)E (2ri ) - V. tan^Y-(2 tan^Y - +ti ) F 
2 o 0 1 X 1 


(6-71) 


where 


L = Log 


Ji 


Tl 


Eo = - W 


-I ^ 

■n X 


dx 


Tl. 

X 
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E^(2ti) = E.(2-n) - E.(2 ti .) 


•n E.(x) - E.(x.) 

F = J dx 


(6-72) 


During the entry phase the altitude variable, t| , is initially 
zero and generally less than 5 when the velocity of the space 
vehicle has been reduced to less than Mach 3. Thus, the exponential 
integrals in Eqs. (6-72) should be evaluated for small t). Expanding 
the exponential term and integrating term by term, gives for E 


/ 

^i 


X 

r; 1 + x+ + • • • 


dx 


X 


= Log (-^) + (ti-ti.) + -rfr (p - r\ ) 
•n. 1 ^ ^ 1 


1 ,3 3, 

+ 37 ^ ) + • • • 


oo 


n n 


= Log(^) + E. 

ri . n= 1 n* n ! 


(6-73) 


F is found by substituting Eq. (6-7 3) into the last of Eqs. (6-72). 

oo 

F = 4- ( Log -3. )^ _ Log E —r~ 

2 ■'li ''li 


oo n n 

E 

” ^ n n! 


(6-74) 


no 



The solution for the moderate and large initial flight path 
angle theory is very accurate and has a wide range of application as 
discussed in the next section. It should be remarked in passing that 
the theory is necessarily limited to large and moderate flight path 
angles simply because the small parameter , 


- 1 

^ ~ — 2 

(3 r V. tan y. 

^ 1 1 

becomes large for small flight path angles. 

Next, the solution will be used to obtain an expression for 
the value of ri at which the maximum deceleration occurs. The 
deceleration in g's is given by 


G 



(6-75) 


The stationary value of G is found from 


dZ ^ 7 dv 

^ ^ d^; 


0 


(6-76) 


With the second of Eqs. (6-53) and the first of Eqs. (6-54), 

v(l - S-q) + e (2 - v) = 0 (6-77) 

To the zero order, e=0, S=S =1, which gives the zero order 

o 

approximation for the altitude at which the peak deceleration occurs 

q = 1 (6-78) 
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The solution can be improved by using the exact equation (6-77) in 
the form 


•n = 1 + e Q(t)) 


where 




, 11 2 
2e tan v. 

1 — 

— —2 " ^ 
(1+ef +e f,) 


0 


Q(ri) = 


- -2 

1 + e + e §2 


(6-79) 


(6-80) 


The transcendental Eq. (6-79) is in a form suitable for the 
application of Lagrange's expansion. Using Eq. (4-45) with r| 
replacing z and Q(p) replacing 4> (p) , the equation for t| becomes 



for small e . To the first order in e the value of r| at the peak 
deceleration is given by 


p = 1 + e 




— . . 2 


(2e - vj tan 


y. + V. 

1 1 


Log (6-82) 


Eq. (6-82) indicates that the altitude at which the peak deceleration 
occurs is sensitive to the initial velocity due to the v^^ Log p^ term 
where p^ is very small. Thus, an improved solution is provided 
which gives the altitude and velocity at the point of greatest drag 

force for any initial velocity and for large and moderate entry angles 

o o 

in the range 5 < - < 90 
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For entry at circular speeds, = 1 and the Log term cancels 
with the expansion of E.(ti.) . 

As suming t) . « 0 , 


oo 

s 

n=l 


n n 
^ -^i 
n* n! 


oo 

E 

n=l 


n 


1 

n* n! 


= 1.31790 


the altitude at peak deceleration can be computed from 


= 


■n smy. 
1 


■ sin Vi 1 + 


1 


1. 31790 + (2e- l)ta 


(3 r tan y. 


V; 


(6-83) 


o o 

which gives very accurate results for 5 < 90 . As a crude 

estimate t) = 1 , in Eq. (6-78), which gives the classical result for 
the critical velocity and the maximum deceleration. 


V 



max 


- 1/2 

e 


0. 60653 


2 


-1 

P r sin y. e 


(6-84) 
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6. 5 Applications of the Ballistic Entry Theories 


For convenience, the three analytic theories for ballistic 
entry of Section 6. 2, 6. 3 and 6.4 will be referred to as the zero 
angle theory, the small angle theory and the large angle theory 
respectively, where the angle is the initial flight path angle. The 
purpose of this section is to discuss the relative merits, ranges of 
applications and accuracies of these theories. 

Fig. (6-1) plots the aerodynamic deceleration versus v for 
zero and small initial flight path angles with pr = 900 using the 
equation 


G 



Z V 


(6-85) 


The dashed line indicates the exact numerical integration of Eqs. 

(6-4) while the solid line represents the zero and small angle theories. 
Inspection of the plot reveals the relative accuracy of the small angle 
theory. In the case of the zero angle theory, however, the solution 
is very accurate, giving 4 digits of accuracy at the maximum 
deceleration of 8. 3 g's. Note that the high accuracy of the zero 
angle theory applies to all values of v from v^^ = 1 down to v = 0. 01 , 
which corresponds to a velocity of about Mach 2. 5 or less. The 
small angle theory is limited in range of application. For larger 
angles, - y. = 10° , there is an obvious discrepancy 
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between the analytic and numerical solutions which disappears as the 

magnitude of the initial flight path angle is reduced. For very small 

values of - , the small angle theory siiffers from the fact that the 

higher order terms begin to diverge. This is due to the appearance 

2 3 

of c / sin Y. in the first order and c / sin y. in the second order 
1 *1 

terms of the analytic solution. Taking just two terms in X from 

Eq. (6-42) and the two largest terms in X from Eq. (6-44) 

2 


e Y, 


e cos Y- 


sin Y. 
1 






2 4 

e cos Y. 

^ 

. 3 

sin Y; 


X^+-l X^l 
12 ^ ^ 36 ^ I 


and then taking the ratio to have a rough convergence test 


i5 

e Y, 


2 

tan Y. 


(X + 3X) 
(3X + 18) 


( 6 - 86 ) 


it is evident that for very small values of - the series is no longer 

convergent. For example, when - X > 1. 385 , or 

V < 0.25 , condition (6-86) is violated. In the figure it is apparent 

that this crude analysis is overly conservative and that actually the 

o o 

small angle theory is limited to the range 2.5 ^ ^ ^ 

In Fig. (6-2) the flight path angle is plotted versus v for the 
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zero and small angle theories. Once again, the zero angle theory is 
very accurate, but the error in the small angle theory is detectable. 
From Eqs. (6-41) and (6-44) , a rough convergence test can be done 
for the flight path angle solution of the small angle theory. 


2 - 


e 


2 

tan y. 


( 4X +9X) 
9(X + 4) 


< 1 


(6-87) 


For - y. < 2 , this condition is violated when X > 3. 254 or when 

1 — 

v< 0. 0386. On the other hand, for - ^ 1° the series diverges 

for X > 0. 954 or when v < 0. 385 . 

The errors due to divergence in the higher order terms cause 
the analytic solution of the flight path angle to be less negative than 
the numerical value. The plot generation of the small angle theory 
is stopped when the absolute error exceeds the absolute error 
obtained if the zero angle theory were used instead. From the figure 
it is apparent that regardless of initial flight path angle, the final 
phase behaves in much the same way as the zero angle entry. 

Fig. (6-3) plots the Log (Z/ Z^) versus v, which by the 


definition of 


Z = 


■ pSC 
2m 


D 


vr = 


pSCj^ r 


2m 




IS 
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= Log - P(r-r^) 

o 


= Q(r - r) (6-88) 

o 

which is the drop in altitude of the entry vehicle. The error in 

both theories is quite small. Numerical computation indicates that 

at maximum deceleration the zero angle theory gives 5 digits of 

2 

accuracy for Log — which amounts to an error of approximately 

o 

3 meters for typical earth reentry. 

In Figs. (6-4) - (6-6) , plots of aerodynamic deceleration in 
g's , flight path angle and Log (Z/ Zo) have been made as a function 
of V for large and moderate flight path angles. The large angle 
theory is very accurate for large angles and has a wide range of 
applications for small angles. The plots show ranges of initial 
flight path angles of -60° down to - 1° to indicate the range of 
applicability in v. The theory actually overlaps the small angle theory, 
providing a definite improvement in the case of y. = - 10° . The 
main defect in the large angle theory is the same as in the others- - 
the higher order terms eventually diverge. The radius of convergence 
of the higher order terms depends on the initial flight path angle, y. , 
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and the initial value of the velocity squared, v^ . This is evident 

from the definition of e = € / v. tan^y. , which must be small to 

1 1 

insure the convergence of the analytic solution of the nonlinear 
differential equations by Poincare's method of small parameters. 

As a result, when - is small, e becomes large, and the range of 
validity of the solution is restricted to small values of the independent 
variable Tj . Furthermore, by the definition q = - 2Z/^p r’ siny. , 


small values of - y. will make q larger at any value of Z. This 

further restricts the application to high altitudes and correspondingly 

to large values of the speed. It is interesting to note from the 

numerical computations that even for very small - y. , for example, 

- y. = 1° , the large angle theory gives very accurate results in the 

highly restricted region, providing 6 digits of accuracy for the 

peak velocity and 3 digits for the minimum flight path angle. An 

error analysis reveals that the magnitude of the difference between 

the analytic solution of v(q) and the numerical integration of the 

2 4 

differential equations (6-4) is approximated by e/ 2 (p r) tan y^ 
during maximum deceleration. For - 60° < Y • £. “ 5 this 
heuristic formula predicts an error which is generally an order 
of magnitude greater than the true error. The effect of the initial 
flight path angle is clearly indicated by the error formula. However, 
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even when the initial flight path angle is -5°, the large angle theory 
gives 3 digits of accuracy at the time of maximum deceleration. 

For = -60° , the theory is extremely accurate, giving 7 digits of 
accuracy for the v(Z). 



^ V 


Fig. 6-1. Variations of G s as Function of v. 

G's = aerodynamic deceleration, v = dimensionless 
velocity squared = V^/ gr and = initial flight 
path angle. The dashed line indicates the exact 
numerical solution while the solid line represents 
the analytic solution from Eqs. (6-85) , (6- 6) , (6-9) , 
(6-19), (6-26), (6-31), (6-34), (6-39), (6-42) and 
(6-43). 
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0.00 0.20 0.40 0.60 0.60 1.00 

V 

Fig. 6-2. Variations of -y as Function of v. 


y - flight path angle, v = dimensionless velocity- 
squared = V^/ gr and y. = initial flight path angle. 
The dashed line indicates the exact numerical 
solution while the solid line represents the analytic 
solution from Eqs. (6-6), (6-9), (6-20), (6-27), (6-31), 
(6-34), (6-39), (6-41) and (6-43). 
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LOG(Z/ZO) 


Fig. 6-3. Variations of Log(Z/ Z^) as Function of v. 

LiOg{Z/ Zq) ~ (Tq-t)/ H = drop in altitude in units 
of scale height 

V = dimensionless velocity squared = V^/ gr and 
= initial flight path angle. The dashed line 
indicates the exact numerical solution while the 
solid line represents the analytic solution from 
Eqs. (6-6), (6-9), (6-19), (6-26). (6-31). (6-34), 
(6-39), (6-42) and (6-43). 


= -15 


V. = -10' 


y. = -5 
1 


Vi=-1‘ 


Variations of G's as Function of v. 

G's = aerodynamic deceleration in g's and v = 
dimensionless velocity squared = V / gr . 

The dashed line indicates the exact numerical 
solution while the solid line represents the 
analytic solution from Eqs. (6-85) and (6-69) 
(6-72). 




y = -6 


Fig. 6-5. Variations of -y as Function of v. 

y = flight path angle, v = dimensionless velocity- 
squared = V^/ gr and y. = initial flight path angle. 
T.he dashed line indicates the exact numerical 
solution while the solid line represents the 
analytic solution from Eqs. (6-69) - (6-72). 







LOG(Z/ZO) 



V 

Figure 6-6. Variations of Log(Z/ Z^) as Function of v. 


Log(Z/ Zq) = (ro-r)/ H = drop in altitude in units 
of scale height, 

V = dimensionless velocity squared = V^/ gr and 
= initial flight path angle. The dashed line 
indicates the exact numerical solution while the 
solid line represents the analytic solution from 
Eqs. ( 6 - 69 ) - ( 6 - 72 ). 
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CHAPTER 7 


CONCLUSIONS AND RECOMMENDATIONS 

A space object traveling through an atmosphere is governed 
by two forces; aerodynamic and gravitational. On this premise 
equations of motion were derived in order to provide a set of univer- 
sal entry equations applicable to all regimes of atmospheric flight 
from orbital motion under the dissipative force of drag through the 
dynamic phase of reentry and finally to the point of contact with the 
planetary surface. The universal entry equations apply to any 
vehicle whether it is a lifting body or a ballistic missile. They 
provide a set of variables which can easily be cast in the form of 
orbital elements or reentry variables so that versatile analytic 
theories can be formulated. This report clearly demonstrates the 
flexibility and wide range of applicability of the universal entry 
equations. 

In the case of orbital motion, a complete uniformly valid 
theory was developed for orbit contraction due to atmospheric drag 
for all elliptic orbits ( 0 < e < 1). A simple model was chosen in 
which the atmospheric density was assximed to vary exponentially 
with altitude, the scale height was assumed constant and the planet 
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and its atmosphere was assumed to be spherical and rotating. 
Rigorous mathematical techniques, such as averaging, Poincare's 
method of small parameters, and Lagrange's expansion, were 
applied in order to obtain a highly accurate, purely analytic theory 
in which the semi-major axis was found as a function of the 
eccentricity. This represents a major breakthrough in the advance- 
ment of analytic theories and should provide ample room for the 
development of more sophisticated theories which include the 
variation of scale height with altitude and the effect of an oblate 
atmosphere and planet. In its present form, the analytic solution 
can be used to develop a computer algorithm which updates the 
scale height at frequent intervals and essentially provides the 
solution for varying p . This method of coupling the theory with 
numerical computation would have the advantage of very short 
computer processing time because the method of averaging has 
removed the need to integrate over every orbit, so that only the 
orbital elements need to be integrated. 

In addition to the orbit contraction, first order solutions 
were obtained for the orientation of the orbit plane, illustrating 
once again the aforementioned techniques. Since the argument of 
the periapsis remains constant and the argument of the ascending 
longitude and the angle of inclination vary very little, the first 
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order solution is adequate for small values of eccentricity. Further 
research work in this area could be done to generate the theory for 
large values of eccentricity and to include the variation of the 
argument of the periapsis. 

In order to give a complete description of the flight of an 
object through a planetary atmosphere, the universal entry equations 
were applied to the problem of ballistic entry. In this case the 
atmospheric density was assumed to vary exponentially with altitude, 
as before, but pr was assumed to be a constant in accordance with 
Chapman's approach and the planet and atmosphere were assumed 
to be spherical and nonrotating so that the entry is planar. The 
great advantages of the universal equations are: they are exact; 
they are free of any restrictive assumptions; and they contain the 
modified Chapman variable Z. This variable permits a single 
trajectory to be solved for specified initial conditions on the velocity 
and the flight path angle that apply to any ballistic vehicle of arbitrary 
mass, area and drag coefficient. With this in mind, the complete 
theory of ballistic entry was developed by analyzing the three typical 
cases: entry from circular orbit, entry from near circular orbit 
and entry at moderate and large initial flight path angles. 

In the final stages of orbit decay, the orbit circularizes and 
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the entry phase is imminent. A first order theory was developed 
which solves for the trajectory of the vehicle during the last orbit. 

The classical Yaroshevskii solution appears as the zero order term 
and is easily derived. It is more laborious to obtain the first order 
term, but its inclusion provides a significant improvement in accuracy. 
With an error in the fourth or fifth digit, it is not necessary to go to 
the second order term in order to improve on the final theory. 

In some orbit decay problems the entry phase may commence 
before complete circularization. In other cases, such as those 
involving recoverable satellites and manned spacecraft, the entry 
phase may be intentionally initiated. For these, the ballistic entry 
theory for small initial flight path angles was developed to solve for 
the modified Chapman Z variable and the flight path angle as functions 
of the velocity. The second order theory provides adequate numerical 
results but it is not sophisticated enough to exhibit all the known 
phenomena which occur during entry from very small initial flight 
path angles. 

There are some cases such as return from the moon, 
atmospheric capture of planetary probes, atmospheric grazing and 
ballistic missile entry in which accurate trajectory analysis would 
be extremely useful. A highly accurate analytic theory for entry 
at large and moderate flight path angles was developed for these 
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purposes which provides the velocity and flight path angle as a 
function of the modified Chapman variable Z, The second order 
theory has a wide range of applicability in the initial flight path 
angle, covering the cases from - Yj > 5° to - < 90° . Any 

initial velocity can be used so that the theory applies to all hyper- 
bolic, parabolic and highly elliptic as well as moderately elliptic 
orbits in the presence of a planetary atmosphere. 

In conclusion, a complete theory has been developed for the 
motion of a ballistic vehicle in a planetary atmosphere. It applies 
to all regimes of flight from the time in orbit through the entry 
phase until the instant of touchdown. The solutions of the governing 
differential equations have been found by mathematically rigorous 
techniques and provide a firm foundation for the development of 
more sophisticated theories. The applicability, versatility and 
generality of the exact universal entry equations have been clearly 
demonstrated. 
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